Edge States, Entanglement Spectra, and Wannier Functions in Haldane's Honeycomb 

Lattice Model and its Bilayer Generalization 



Zhoushen Huang and Daniel P. Arovas 

Department of Physics, University of California, San Diego 
(Dated: March 5, 2013) 

We study Haldane's honeycomb lattice model and a bilayer generalization thereof from the per- 
spective of edge states, entanglement spectra, and Wannier function behavior. For the monolayer 
model, we obtain the zigzag edge states analytically, and identify the edge state crossing point kc 
with where the / = | entanglement occupancy and the half-odd-integer Wannier centers occur. A 
continuous interpolation between the entanglement states and the Wannier states is introduced. We 
then construct a bilayer model by Bernal stacking two monolayers coupled by interlayer hopping. 
We analyze a particular limit of this model where an extended symmetry, related to inversion, is 
present. The band topology then reveals a break-down of the correspondence between edge and 
entanglement spectrum, which is analyzed in detail, and compared with the inversion-symmetric Z2 
topological insulators which show a similar phenomenon. 



1. INTRODUCTION 

Haldane's honeycomb lattice model^ has provided a 
fertile paradigm for topological band structures in the 
absence of net magnetic flux. Prior to Haldane's work, 
Thouless et al? (TKNN) showed how a tight binding 
model with uniform rational flux per plaquette, i.e. the 
Hofstadter model, yields a topological band structure in 
which each energy band n is classified by an integer topo- 
logical invariant C„, known in mathematical parlance as 
a Chern number. (In the continuum limit, where the 
flux per plaquette is infinitesimal, the TKNN bands be- 
come dispersionless Landau levels.) The essence of the 
Haldane model lies in its inclusion of complex second- 
neighbor hopping amplitudes, so that the model breaks 
time-reversal (T) symmetry even though the net mag- 
netic flux per plaquette vanishes, which allows for the 
existence of topological phase with band Chern indices 
±1. 

One of their hallmarks of topological phases is the ex- 
istence of gapless edge modes interpolating the bulk gap 
in the presence of an open boundary. The number of 
such edge spectral flows as functions of the momentum 
parallel to the boundary is the same as the total Chern 
index of the bands below the gap, as first elucidated by 
HatsugaP. Kane and Mel^ later generalized Haldane's 
model by introducing spin and treating the (now purely 
imaginary) second neighbor hopping as a spin-orbit cou- 
pling. T-preserving perturbations are also allowed. The 
Kane-Mele model is T-invariant and cannot have a quan- 
tum Hall effect. The bulk topological property is instead 
described by a Z2 topological indejiP^. Remarkably, at 
half filling, while the total Chern index is zero, the gap- 
less edge spectrum persists due to time-reversal symme- 
try. A topologically trivial band insulator, by contrast, 
would have no edge spectral flow at all. 

The gapless edge spectral flow of topological insulators 
is one of the real space manifestations of their bulk topol- 
ogy. A similar spectral flow can be observed in the quan- 



tum entanglement spectrum of the many-body reduced 
density matrix obtained by partitioning the system along 
a translationally-invariant boundary^. For noninteract- 
ing fermions, the spectrum of the reduced density matrix 
itself corresponds to that of a noninteracting 'entangle- 
ment Hamiltonian' determined by the one-body correla- 
tion matrix of the original system® There are however 
exceptions to the entanglement and edge spectra corre- 
spondence. For example, the entanglement spectrum has 
protected midgap modes for systems with inver sion (X) 
symmetry even if the edge spectrum is gappecP^^ by 
6.17. breaking the T symmetry in quantum spin Hall ef- 
fect (QSH) systems. In certain cases, one also has to 
tune the boundary conditions for a system with nontriv- 
ial topology in order for its energy edge modes to be 
gapless^^, while such tuning is not required to observe 
the entanglement spectral flow. We shall see similar dif- 
ferences in our study of the Haldane models. Thus in cer- 
tain sense, the entanglement spectrum reveals the bulk 
topology better than the Hamiltonian's edge spectrum. 

The band topology can also be considered from a 
Wannier functior^ point of view. While in higher 
dimensions, the construction of exponentially localized 
Wannier functions is precluded for band insulators with 
nonzero Chern number^i^, the system is effectively one- 
dimensional if one specifies the momentum along the 
edge/surface, for which the Wannier states are well de- 
fined. The Wannier functions are then localized along 
strips or planes parallel to the ed ge. Several recent stud- 
ies of topological insulatorJlSlIIl have invoked the Wan- 
nier states in their analysis. The Wannier centers are 
shown to exhibit a spectral flow similar to that of the 
entanglement spectrum, and the topological information 
can be visually extracted from their flow pattern. Mathe- 
matically, the deviation of Wannier centers from the cor- 
responding unit cells are eigenvalues of the Wilson loop 
operator VV^®. It is interesting that VV is an object de- 
rived purely from the bulk (for translationally invariant 
systems) and is hence faster to compute, yet its eigen- 
values have a real space interpretation similar to the en- 
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tanglement spectrum: when the Wannier centers migrate 
from one unit cell to its neighbor, there is a correspond- 
ing flow in the entanglement spectrum if the particular 
unit cell boundary is used as the entanglement cut. The 
entanglement spectrum is thus a coarse graining of the 
Wannier centers with an emphasis on the real space be- 
havior near the entanglement cuiP^. 

In this paper, we study the Haldane honeycomb lattice 
model and a bilayer generalization thereof, both from a 
real space perspective, combining the analysis of Hamil- 
tonian edge states, entanglement spectra and Wannier 
center flow. We flrst present an analytical solution of 
the monolayer zigzag edge modes, identifying the k point 
where the the two edge modes cross. This plays the role 
of one of the T-invariant k points in the QSH models. We 
flnd that at the same k point, there is an entanglement 
occupancy mode flxed at / = ^, and the correspond- 
ing Wannier centers reside exactly in the middle of two 
neighboring unit cells. We show that a common origin 
underlies this coincidence. We then extend to a bilayer 
model by Bernal-stacking two monolayers with vertical 
interlayer hopping, as in bilayer graphene. With a partic- 
ular parameter choice, the bilayer model exhibits some- 
thing similar to the Z2 inversion-symmetric topological 
insulators (ITI) in that the edge spectrum is gapped yet 
the entanglement spectrum has protected f — \ modes. 
However, it cannot be explained in the ITI framework 
due to the lack of I symmetry (the protected modes do 
not occur at the inversion-symmetric k points). We show 
that the nontrivial topology is a consequence of a related 
T* symmetry to be detailed in the text, and derive an ex- 
pression to compute the topological index, which is the 
winding number of one branch of the Wannier centers. 
We further conflrm this numerically by adding in T*- 
preserving perturbations to both the bilayer model and 
an ITI model studied in Ref. [121 



2. MONOLAYER HALDANE MODEL 

We first briefly review the Haldane model, which is a 
tight binding model on the honeycomb lattice, described 
by the Hamiltonian 



n 



(1) 



The hopping amplitudes t^j are nonzero only for nearest 
neighbor (NN) and next-nearest neighbor (NNN) hop- 
ping. For NN hops, tij = t is real. For NNN hops, 
tij = tge^'^'^, where s = 1 if the hopping is parallel 
to fli, s = 2 if parallel to (22, and s = 3 if parallel to 
aa = a2 — ai . The sign of the phase is according to the 
arrows in Fig. [TJ and is taken to be positive for clockwise 
hops within each hexagonal unit cell. In Haldane's orig- 
inal model, ti — t2 = t^. Setting these amplitudes to be 
different breaks the three-fold rotational symmetry. The 
Semenoff mass is m and —m for A and B sublattices 
respectively, which breaks inversion. 




FIG. 1: (Color online) Haldane's model. A (red) and B 
(white) label sublattice sites, and the boxed pair represents a 
unit cell. Primitive lattice vectors are chosen to be ai and 02 
as shown. Second neighbor hopping between same sublattice 
sites picks up a phase of along the arrowed directions and 
— (j) in opposite directions. Horizontal box encloses a zigzag 



In the bulk, the Fourier transformed Hamiltonian is 

H{k) ^uj{k) + B{k)-a, (2) 

— (crxjCyjaz) are the Pauli matrices in the "isospin" 
degree of freedom, where A and B are isospin up and 
down, respectively, and 

Bx = —1 — cos fci — cos k2 (3) 
By = — sin ki — sin k2 

Bz = m + 2 sin 0[ti sin fci — t2 sin k2 + t^ sin(fc2 — fci)] 

UJ = —2 COS(f>[ti cos fci + <2 COS fc2 + ^3 C0s(fc2 — fci)] 

Here, fcj — k ■ Hi e [0, 27r] are the Bloch phases along 
the two primitive vectors. The bulk topology is char- 
acterized by the Chern number C of the upper band, 
which is the winding number of the unit vector B{k) 
over the Brillouin zone. That is to say, if by varying 
(fci,fc2) over the first Brillouin zone, 13 = B/\B\ covers 
the unit sphere once (|C| = 1), then the system is in 
its topological phase, otherwise it is in its trivial phase. 
Equivalently, in the topological phase the origin is inside 
the surface swept out by B, while in the nontopological 
phase it lies outside. The topological phase transition 
thus takes place when the origin is crossed by B at some 
fc points, where the gap A = 2\B\ will vanish. Following 
eqn.[3j this can only happen at the graphene Dirac points, 
(fci,fc2) = ±(27r/3,-27r/3) = K±, where B^ ^ By ^ 0. 
The corresponding B^ values necessarily have opposite 
signs in the topological phase (so that the origin is en- 
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closed), i.e!^, 

H-V3|(t,+t, + .3)sin0||<O 'X!^'''- (4) 



In the trivial phase, the eigenstates just below the gap at 
both Dirac points are entirely concentrated on the same 
sublattice. In the topological phase, however, they are 
concentrated on opposite sublattices. 



2.1. Zigzag edge states 

One way of solving the edge spectrum of tight binding 
models is to use the transfer matrix formalism, following 
Hatsugai's investigatior!^ of the Hofstadter problem^°. It 
is worthwhile to first think about its application in the 
Haldane model. In the Hofstadter problem, the system 
is immersed in a uniform magnetic field with rational 
fiux 2TTp/q per (square, say) plaquette, and the lattice 
vector potential is periodic on the scale of the magnetic 
unit cell, which consists of q structural cells. As a re- 
sult, a giV-step transfer matrix A4qN can be broken up 
into a product of N identical matrices each equal to a q- 
step transfer matrix, MqN = {Mq)^ , thus the solutions 
oi Mq comprise a special set of the solutions of MqN- 
Physically this means the edge spectrum of a system with 
qN — 1 structural cells in the open direction is identical 
to the full spectrum of a system with q — 1 structural 
cells. Thus numerically, one only needs to diagonalize a 
(g — 1) X (g— 1) Hamiltonian to find the edge spectrum of 
a (qN — l) X {qN— 1) Hamiltanion. Clearly the method is 
most efficient in a situation where the transfer matrix has 
such a decomposition, e.g., a graphene sheet in magnetic 
ficld^^. The Haldane model has no macroscopic magnetic 
field (an essential virtue of the model) , thus the transfer 
matrix formalism yields little advantage over directly di- 
agonalizing the full Hamiltonian. Still, one may employ 
it to analyze the Riemann sheet structure of the complex 
energy, which is studied in Ref. 22, but that is not our 
focus here. 

A useful feature of Hatsugai's solution is that it can be 
written as a direct product of an A^-component real space 
part, corresponding to the magnetic unit cell coordinate, 
and a g-component internal space part, corresponding to 
the lattice points within each magnetic unit cell. (See, 
e.g., the appendix of Ref. This is by no means a 

general form for edge states. All Bloch states on the 
other hand have such a decomposition. What we found 
for the Haldane model is that the edge states in the case 
of a zigzag edge can also be direct-product-decomposed, 
with the following caveats. First, while the real space 
part in the Hofstadter model has exact exponential de- 
pendence on the coordinate (equivalent to an imaginary 
Bloch wavevector), which is due to decomposition of the 
transfer matrix, this is not so in the Haldane model (nor 
is this surprising since the macroscopic magnetic field is 
zero). Second, the boundary condition used in the Hof- 



stadter model corresponds to an open edge. In the Hal- 
dane model, the boundary condition must be tuned self- 
consistently to conform with the direct product Ansatz. 
Only those boundary conditions with vanishing magni- 
tude will correspond to an edge state at an open bound- 
ary, as opposed to, say, an enhanced-tunneling boundary. 
We now proceed to solve for the zigzag edge states. 



2.1.1. Twisted-boundary Hamiltonian and gauge 
transformation 



The zigzag edge is parallel to di (see the horizontal box 
in Fig. [ij, hence ki — k ■ di is a good quantum number. 
Assume there are N unit cells in the a2 direction. At each 
ki , the effective 1-D system is described by the following 
Hamiltonian, 

N 



n,n' — l 



(5) 



Here a\ and h\ are creation operators on the A and B site 
of the n**^ unit cell, respectively. The coefficient matrices 
ii'^*^ connect sites on the same sublattice, and R connects 
different sublattices. Their nonzero matrix elements are 
given by 



V* hi V,, 



R = 



\Vi 

(VI 

P2 Pi 



V* hi V, 
V* hj 



(6) 



P2 Pi 

P2 Pi 



(7) 



The individual matrix elements can be obtained from 
Fourier transforming the bulk Hamiltonian eqn. [2j 

'hi pi 
J>i 



J2H{ki,k2 



P2 V2 



(8) 
(9) 



with 



hi = m — 2ti cos((/) 
/i2 = — rn — 2ti cos( 



fci) , 
-fci) , 



vi = -t2e-"'' - tse^t'e-'''^ 



V2 = ~t2e 
Pi = - 



-iki 



P2 



(10) 
(11) 

(12) 
(13) 
(14) 
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We note that swapping subscripts 1 and 2 on the left hand 
sides yields a Hamiltonian with the so-called bearded 
edge. The method we describe below applies to both 
types of edge. 

The following gauge transformation makes both pi and 
P2 real, 



by which Vi — >■ Vi e*'^i/^ and 



pi — pi e'*"'^/^, viz.. 



vi -{t2 + h) cos{(t> - + i{t2 ~ h) sin( 



(15) 



2 



V2 



-{t2 + h) cos{(t> + ^) - i(i2 - ts) sin((/) + 



pi —T' -2 COS ^. 



(16) 



In eqns. [6] and [?[ a twisted boundary conditiorP^ is 
used: 



V = p Uv = 



Vl Vi2 
V21 V2 



(17) 



p is a real number controlling the "tunnelling strength" 
between the two edges, and ?7 is a unitary 2x2 ma- 
trix that describes an "isospin-dependent" phase twisting 
over the boundary. For an open boundary, /? — )• 0. For 
periodic boundary conditions, p = 1 w ith U = I (without 
the gauge transformation of eqn. 15) or [/ ~ e^'^'^'^/^I 
(with eqn. 15), where I is the 2x2 identity matrix. In- 



troducing twisted boundary condition may seem to over- 
complicate the situation, but as we shall see it allows us 
to make progress toward an analytical solution. 
The eigenvalue problem can now be written as 



K<''^4'a) +R\M =£\M , 



(18) 



where \iPa) and \iPb) are the "wavefunctions" of the A 
and B sublattices, both of which are TV-dimensional col- 
umn vectors. 



2.1.2. Edge state Ansatz and energy 

We look for solutions of the form \iIja) = and 
I ■0b) — ■ 111 terms of the direct-product decomposi- 
tion discussed earlier, \^) is the real space part and {\) 



is the internal space part. Eqn. 18 now becomes 



(if (1) - 
(Aif(2) 



\R 



-em . 

A£)|^> 











(19) 
(20) 



A sufficient condition for both equations to be satisfied 
is that the coefficient matrices are proportional clement 
by element, 



hi + Xpi — £ 



«1 + 



A ipi — e V2 + X ^P2 



r. (21) 



This gives, at each value of fci, two equations (the ratio 
r itself being yet undetermined) for the two unknowns A 
and e. The solutions are 



A. 



V1V2 - V2 



vI-pI± ^(P-A\viV2\ 



(22) 



_ d±^cP-A\v,V2\^ 
r± = r(A±) = — — ^ , (23) 



2\V2\^ 

r±{p2h2 - 2piRet;2) - (^2^1 - IpiRevi) 
P2{r± - 1) 



(24) 



where Re indicates the real part, ± denotes the branch 
of solution, and 



+ V2V1 - P2 & 



(25) 



See Appendix [Xjfor details regarding derivation. 

Eqn. [23] becomes singular when V2 = 0, which could 
happen for the zigzag edge if ^2 = is and fci = tt — 2(f>. 
For this particular parameter set, one can readily verify, 
by Taylor expansion in V2, that 



X+ =vi , r+ = -v-i , e+ 
Al^ = W2 ^ , rZ^ "2 



h2vl + 2pivi + hi 
1+^ 



1-2-^0, £_ ft.2 . (26) 



These results also hold both for graphene (m = 0) and 
for boron nitride (to 7^ 0). In both cases, second neighbor 
hoppings are turned off, rendering vi = V2 = 0. Clearly, 
A+ = r+ = Al"^ = rZ^ = 0, and e± — ±m. As will be 
shown in Appendix [A] p — if |fci| > 27r/3, so solutions 
there correspond to edge modes with open boundary. 

A natural question arises regarding the reality of the 
energy eqn. 24 As long as we can find wavefunctions 
complying with the Ansatz lips) — ^{iPa), e± will be 
eigenvalues of a Hermitian matrix, and hence real. This 

) I . Thus by 
yields real solutions, with 



implies that r± are also real (c/. eqn 
I Ansatz 



23 



eqn. 

choice of p [/, provided the discriminant satisfies 
A = -4|wiU2p > . 



some 



(27) 



Note that this condition is valid for all ki. 

Although real solutions exist for all ki with A(fci) > 0, 
they do not necessarily correspond to open boundaries. 
Normally, wavefunctions are solved after fixing boundary 
conditions (p and U), but here, we enforced a particular 
form of solution, which will not be consistent with arbi- 
trary p U. Instead, the matrix p J7 is to be determined 
self-consistently from the Ansatz, and is in general ki- 
dependent. This is discussed in detail in Appendix 
For now, we just note that only when p — > will the 
solution be valid for open boundary. Clearly, as p varies 
with ki, the transition from an open boundary solution 
to that of an "enhanced tunnelling" boundary will hap- 
pen when p = 1, at which point the edge solution merges 
into the bulk (|'0) becomes extended instead of localized) . 

Fig. [2] shows the Ansatz solutions (colored curves) 
comparing with the open boundary spectrum (circular 
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ME+) 

Re(B_ 
Open boundary 

I \ 

-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

ki/iT 

(a) Ansatz vs. open boundary spectrum 




(b) auxiliary quantities 



FIG. 2: (Color online) Zigzag edge modes for [ti,t2,ts] = [0.3,0.4,0.5], m = 1.4, = O.Svr, A*' = 40. (a) Circular dots are 
obtained from diagonalizing an open boundary Hamiltonian, eqn. [Sjwith p = 0. Colored curves are Ansatz solutions, eqn. |24| 
Gray vertical lines mark the edge level crossing points given by eqn. |30[ (b) Auxiliary quantities. A is the discriminant in 
eqn. |27| only A > yields real solutions of edge energy e± and is physical. r+: the ratio of eqn. [21] corresponding to the 
is real if A 



£+ branch. It is real if A > 0. r_ can be obtained from eqn 

strength as defined in eqn. [5] These are solved in a self-consist fashion in Append 
"enhanced-tunnelling" boundary. The transition p — 1 corresponds to bulk modes 



A3 and is not plotted here. 

P- 



p{e±): the inter-edge tunnelling 
> for open boundary, p > 1 for 



dots) in |2(a)[ and auxiliary variables A, r{e+) and p{e±) 



2(b) 



Parameters are chosen to exhibit most of the 
possible scenarios, [^1,^2,^3] — [0.3,0.4,0.5], m — 1.4, 
(j) = O.Stt, iV = 40: 

(1) A < for |fci| < 0.125. In this region the Ansatz 
does not yield real energy solutions. 

(2) For |fci| > 0.125, the Ansatz yields physical solu- 
tions. One then computes p self-consistently; p I 
means open boundary, whereas p > 1 means enhanced- 
tunneling boundary. The transition happens when the 
open boundary edge modes merge with the bulk bands. 
The e+ branch (red curve) has open boundary for fci e 
[-TT, -0.627r] U [-0.177r,-0.147r] U [0.627r,7r], while the 
e_ branch (blue curve) has open boundary for ki G 
[— TT, — 0.487r] U [0.497r, tt]. In these regions, the Ansatz so- 
lutions overlap with the open-boundary numerics (filled 
circles). Note that the £+ branch briefly becomes open 
boundary in [— O.ITtt, — 0.147r]. Without the Ansatz so- 
lution, one would have taken it to be part of the bulk 
spectrum. 

(3) Within the physical regime (A > 0), the two branches 
cross twice, marked by the two vertical gray lines, one 
with p{e±) — >■ and the other with p(e±) > 1. In both 
cases p{s+) = p{s-). These two edge crossing ki points 



are described by eqn. 30 which will be discussed in the 



next section. While only the one with p — is relevant 
for the open-boundary edge spectrum, we shall see in the 
next section that both have geometrical significance and 
will be reflected in the entanglement spectrum and Berry 
phase flow. 



2.2. Topological signatures in edge, entanglement 
and Wannier spectra 

A gapless edge spectral flow is one of the most con- 
spicuous real-space manifestation of nontrivial topology 
of band insulators. But as discussed in the introduction, 
it sometimes fails to reveal every topological difference a 
band insulator can have from its atomic limit^^, which 
the entanglement spectrum can capture. In the non- 
interacting fermion case, the entanglement occupancy 
spectrum is the eigenvalues of a sub matrix G of the one- 
body ground state projector ^^ESl^ with the dimension 
of G given by the entanglement cut. The eigenvalues of 
Q itself are just Os and Is, but for a system with non- 
trivial topology, the eigenvalues of G exhibit a spectral 
flow from to 1, as a function of the momentum kj_ 
along the entanglement cut. The reason why the en- 
tanglement cut would induce such a flow can be under- 
stood intuitively in terms of the Wannier states. One 
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V 


Lower band 




Upper band 



(a) Periodic-boundary entanglement 
occupancy 



(b) Periodic-boundary entanglement 
quasi-cncrgy 



-1 -0.5 0.5 1 

Ai/ir 

(c) Deviation of Wannier centers from unit 
cell boundary 



FIG. 3: (Color online) Coincidence of level crossings in entanglement and Berry phase flows, (a): Entanglement occupancy, (b); 
Entanglement quasienergy. (c): Deviations of Wannier centers from their corresponding unit cell boundary; mathematically this 
is the Berry phase; plotted for both the occupied and unoccupied bands. Periodic boundary conditions are used. Parameters 
used here are the same as in Fig. [5] [ti,t2,i3] = [0.3,0.4,0.5], m — 1.4, cf> — O.Stt. Entanglement spectra are computed for the 
lower half system with M — N/2 — 20 unit cells in the a2 direction, where A*' = 40 is the number of unit cells of the full system, 
and Fermi energy is placed in the gap, — 0.6. Berry phases are computed using 100 steps in the k2 integration. Notice the 
similarity between the entanglement occupancy and the Berry phases. This is due to the former being a coarse-grained version 
of the latter in the sense of Ref. 1191 Notice that in all three plots, level crossings occur at the same kc points given by eqn. |30| 
Although the open boundary edge modes only cross at one of the two kc, the entanglement spectra and the Berry phases cross 
at both kc, but with different /, e and 7 values, f — e — and 7 = tt at the open-boundary edge crossing point. At the 
other kc point, 7 = 0. Notice also that in the quasienergy plot, level crossings are not restricted to the central two levels but 
extend all the way to big I7I values, and are pinned at the same kc values. 



can always recombine the Bloch states that constitute Q 
(assuming the system is periodic) into a set of spatiaUy 
localized Wannier states. If a Wannier state resides in 
either half of the partition, it has almost perfect projec- 
tion onto that half, and the corresponding entanglement 
occupancy is very near or 1. If on the other hand the 
entanglement cut passes right through a Wannier state 
(the exact meaning of which we shall discuss in ^ 2.2.3 1, 
then it has significant projection onto both halves, cor- 
respondingly the entanglement occupancy is near i. In 
the topological phase, the Wannier states themselves flow 
with respect to fc ^ i^S nis ^ ^j^^g when one such state passes 
through the entanglement cut, a corresponding entangle- 
ment flow arises. One can see that the position of the 
Wannier states, or the Wannier spectrum, is closely re- 
lated to the entanglement spectrum, and that a nontriv- 
ial topology underlies the f = \ entanglement occupancy 
mode, and a half-odd-integer Wannier center. It is thus 
interesting to note that in the Haldane model, the edge 
crossing point, the half occupancy mode and the half- 
odd-integer Wannier center all coincide at the same fci 
point. See Fig. |3] In this section we shall study the rea- 
son underlying this coi nciden ce. We mention that there 
are several other work^^SEll studying the charge polar- 
ization of the Haldane model, from the Wannier function 
perspective. 



2.2.1. Edge modes crossing points 



With the edge solution in j ]2.1[ we can identify the 
exact location of these crossing points. The condition for 
£4- = e_ is that 



P2{hi 



(28) 



see Appendix \K\ for derivation. For the zigzag edge, this 
implies 



m + 2 sin(/)sin kciti -f ^2 + ta) = 



(29) 



whence 



m 



2{ti+t2+h)sin(l) 



(30) 



TT + sm 



2{ti+t2+h)ii\TL(j) 



Here denotes the values of fci where the two edge 
modes are energetically degenerate. The bearded edge 
solution is obtained by switching the 1 and 2 suffixes in 



Pi, hi and Vi, for which eqn. 28 implies 



(m -|- 2ti sin (j) sin fee) cos ■ 



+ {t2 + h)sin<f)sm^ =Q. (31) 
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This could be recast as a cubic equation for tan(fcc/2), 
but this does not afford a particularly simple closed form 
solution. 

One can see from Fig. [2] that only one of the two fee 
has p — ^ and corresponds to an open boundary edge 
mode. The other one is of an enhanced-tunnelling bound- 
ary; in some parameter settings it even lies in the region 
A < where the Ansatz solution is complex. However, 
the enhanced-tunnelling kc is still special in the entangle- 
ment and Wannier center flows, as can be seen in Fig. |3] 
What is the significance about these edge crossing points? 
Recall that the bulk Hamiltonian H{ki, /C2) maps each k 
to a B{k) vector. Fixing ki while varying fc2 will drive 
the B vector along a closed curve in 3D. It turns out that 
at both kc, this curve lies on a plane passing through the 
origin. To see this, we note that at the edge crossing 
point, the 2x2 blocks h{ki) and v{ki) are related via 



P2 



where e, 

Hamiltonian at the edge crossing points is then 



e+(fcc) = e^{kc) (cf. eqn. A8). The bulk 



H{kc,k2) 



h 



Pi 

P2 



Ak2 



H.c. 



(33) 



and the corresponding B{kc, k2) is 



i?2 = ( — -f cosfc2 I Re (wi - V2) - sinfc2 Im (ui - V2) , 
\P2 J 

Bx = Pi + P2 COS k2 , By = p2 sin ^2 , (34) 

note in particular that all components are independent 
of TO, which is the term that breaks the inversion of the 
two sublattices. It is then easy to check that 

B(fcc, fc2) = (i - Re(g) z) + By{y~ \m(q) z) (35) 
where 

V2 - Vi 



q = 



P2 



(36) 



is independent of fc2. The path of _B(fcc, fc2) is thus copla- 
nar and normal to the vector 



Re{q) X + \m{q) y + z. 



(37) 



An unrestricted {Bx,By) pair can describe any point on 
the plane. Clearly the origin itself is on the plane. The 
actual path of B{kc, k2) is restricted to those allowed by 
eqn. [M] 

An interesting observation is that in graphene, the bulk 
2x2 Hamiltonian is always off-diagonal, thus it is copla- 
nar at any ki value. The origin is inside the path of B 
for I fell > 27r/3, and outside otherwise, hence the well 
known result that its two zigzag edge modes are degen- 
erate at e = for > 27r/3. For the bearded edge, the 
degenerate edge modes connect the two Dirac points in 
the other way, namely with |fci| < 27r/3. 



2.2.2. Integer and half-odd-integer Wannier centers 

For a general one-dimensional periodic system, or 
higher-dimensional system with k^ specified, the Wan- 
nier centers can be defined as the nonvanishing eigenval- 
ues of the band-projected real-space operator 



Q{k^)Rg{k^) 



(38) 



where ^/(fc_|_) is the filled band projector (or sum of pro- 
jectors for multiple bands) at fc_L, R = Y ®1 in which 
Y = diag(l, 2, . . . , N) measures real space (i.e., unit cell) 
coordinates, N being the number of unit cells in the lon- 
gitudinal direction, and I is the g x g unity acting on the 
q-dimensional internal space. The corresponding eigen- 
states are defined as the Wannier states^. The mono- 
layer Haldane model has only one band occupied at half 
filling. When Q contains only one band, the eigenvalues 
areP 



(32) of eqn. 38 



A/(fc^) = ^+/ , /=l,2,...iV (39) 

where / labels the unit cell, and 7(fcj^) is the Berry phase 
of the band, 

-l{k^)= (dk.Mk..) , A{k..) = i{tp\d, . (40) 



The /cj_ dependence of A is suppressed. A is the fcy com- 
ponent of the Berry connection vector at fixed k^ . Thus 
the offset of the Wannier centers from the unit cell bound- 
aries are uniform and are given by the Berry phas^^SI 



In the Haldane model, k 
lower band projector is 



ki and fc,. 



g(fci)==^|vl/(fc))(vl/(fc)| 



|*(fc)) = |fc2)®|^(fcl,fc2)) 



The 

(41) 
(42) 



where the internal part |i/'(fci,fc2)) is the lower band 
eigenstate of the bulk Hamiltonian eqn. [2] and the real- 
space part \k2) is the Bloch wave in the 0,2 direction, 
{y\k2) — e^^^y I^/N. The upper band projector is simi- 



larly defined. We plot the Berry phase in Fig. [3(c) where 
the red and blue curve correspond to the lower and upper 
band, respectively. We will show below that coplanarity 
of B at kc fixes the Berry phases there to be either tt or 
0, depending on whether or not the origin lies inside the 
path of B. In fact, this holds for any two-band model 
with coplanar k points. Recall H{k) — cl>(A:)I + B{k) ■ a. 
Let [B{k),d{k),(p{k)] be the spherical coordinate of 



(43) 



B{k), then the lower band eigenstate is 



sm 2 
cos ■ 
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whence 



and 



A{k2) = -i(l + cosi?) 



dip 



(44) 



7(fci) = -\ dip {1 + COS'S) 



(45) 



v(o) 



where i^(0) = ip{ki,k2 = 0) and (/3(27r) = ip{ki,k2 — 27r). 
At the edge crossing fci, i3 is comprised of points on a 
great circle (since the origin is on the plane), so one can 
always rotate the internal space to a frame where the 
plane normal n (eqn. 37) coincides with the z axis, then 
cos ?? = everywhere on the coplanar path which is now 
lying in the new xy plane. The Berry phases at coplanar 
ki points are simply 27r times half the (negative) winding 
number w of the path of -B, 

7(fcc) = -w , (46) 

Note that this rotation amounts to applying the follow- 
ing unitary transformation to the band projector Q in 
eqn.|351 

g^g^ugu^ , u^i(E)D(n) (47) 

where D{n) is the aforementioned internal space rotation 
that depends only on n. Since URW — R, both ORG 
and ORG have the same spectrum, so the resulting Wan- 
nier centers and Berry phases are independent of whether 
Q or Q is used. 

For the Haldanc model, the topological phase has 
Chern number \C\ = 1, so the winding number at any 
ki is at most |w| = 1, hence "f{kc) is either ±7r (ori- 
gin inside) or (origin outside), which is what we see in 
Fig. |3(c) In the non-topological phase, 7 at both fee will 
be 0. 



2.2.3. Entanglement Half Occupancy Mode 

As discussed in the introduction, the entanglement 
spectrum can be associated with a coarse-graining of the 
Wannier centers. By coarse-graining, we mean replacing 
the real-space operator R — Y ^1 with the projector 



fiy,M) 



Pm = fiY, M) ® 1 , 

\ iiy <M 
if y > M 



(48) 
(49) 



where M is an integer. In so doing, all Wannier center 
flows except those between unit cells M and M + 1 are 
suppressed. Since Q and Pm are both projectors, QPmQ 
and PmQPm share the same spectrum (cf. Appendix [Pj). 



The latter is nothing but the restricted correlation ma- 
trix whose eigenvalues are the entanglement occupancy 
spectrum, and the integer M corresponds to the entan- 
glement cut. The entanglement eigenstates, j^*), are pro- 
jections of eigenstates \^) of QPmG onto the half space 
y<Mhy Pm. i-e. |*) =Pm\^)- 

We found that at the edge crossing point fcc, in the case 
of odd winding number, the periodic boundary entan- 
glement occupa ncy sp ectrum has two modes intersecting 



icy sp c 

at / = ^ (Fig- 3(a)|, which corresponds to a zero en- 
tanglement quasi-energy (Fig. [3(b)] ). It is evident from 
Fig. |3(b)| that the entanglement quasi-energy spectrum 
has a particle- hole symmetry for all ki. This is a conse- 
quence of the periodic boundary condition: there are two 
independent flows in Figs. 3(a) and |3(b)] - one upward, 
one downward - because by using periodic boundary con- 
dition, an entanglement cut creates two new boundaries 
in each of the half systems. These two flows intersect at 
both fcc points. Had we started with a Q for an open 
boundary system (in the 02 direction), then the entan- 
glement cut would only create one new boundary to each 
of the half systems, and as a result, in each half sys- 
tem, only one of the two flows, which corresponds to the 
edge created by the entanglement cut, would survive. In 
that case, there would be no particle-hole symmetry in 
the entanglement quasi-energy spectrum for arbitrary fci , 
but only at the edge crossing kc points, where the entan- 
glement spectrum still exhibits particle-hole symmetry, 
since switching the 0,2 boundary condition from periodic 
to open merely removes the double degeneracy at these 
points. 

In fact, one can argue that in an open boundary sys- 
tem, the f = \ level is a consequence of this survival of 
entanglement particle-hole symmetry at the edge crossing 
point. Assume the number of unit cells in the full system 
is N and N is even. In an open boundary system, there 
are two edge levels, thus each bulk band has — 1 bulk 
levels. At the edge crossing point, both edge levels are 
either below or above the Fermi energy, so the rank of 
Q is A^ ± 1 which is always odd. We may always pick 
the one half of the bipartite whose size M > rank(C/), so 
that rank(G) = rank(CJ). Thus with the entanglement 
particle- hole symmetry at the kc point, there must be an 
odd number of entanglement levels fixed at / = ^. 

With a periodic boundary, the half occupancy modes 
will appear in the thermodynamic limit, and for small N, 
there is avoided crossing between the two flows. How- 
ever, the degeneracy is exact for arbitrary even N and 
M = N/2 when t2 = t^. We prove this in Appendix [b| 
Below, we wish to discuss it following the coarse-graining 
perspective. 

The effect of inserting R and Pm in between two band 
projectors G is to resolve, in real space, the Bloch states 
that constitute G- As a result, both the Wannier states 
and the entanglement states are spatially localized. In 
this respect, there exists a family of such real-space re- 
solvers Rp, parameterized by an "inverse temperature" 
/3, that smoothly interpolates between the Wannier and 
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entanglement limits, 

tanh 



i/3(f-i(7V + l)) 



Rp = \ 



tanh [|/3(A^ - 1)] 
{N -l)fp(Y) + N + 1 



(50) 
(51) 



Here, fp is obtained by rescaling and shifting the Fermi 
distribution. The extremal cases are 



R 



tS-s-O 



Y 



N- 1 



P. 



N/2 



(52) 
(53) 



where I — P/v/2 is a real space projector onto the y > N/2 
half. Thus by varying /3, the spectrum {r^} of GRjs Q 
morphs from the Wannier center spectrum to the entan- 
glement spectrum, providing a one-to-one mapping be- 
tween the Wannier states and the entanglement states. 

At the coplanar point with winding number w = 1, 
two ra levels are pinned at {N + l)/2 regardless of the 
value of /3. This is shown in Fig. |4] The /3-independence 
suggests a unifying physical origin underlying these lev- 
els. We have shown before that, mathematically, the rea- 
son for the Wannier spectrum (/? — 0) to exhibit such 
levels is because the Berry phase of the occupied band, 
which corresponds to the uniform deviation of all Wan- 
nier centers from the integers labeling their unit cells, is 
precisely tt at the w = 1 coplanar point. It is instructive 
to reflect on the physical implication of this result. For 
convenience, let us rotate to the internal basis where the 
path of B lies in the xy plane, cf. discussion leading to 
The Berry connection and its cumulation along 



46 



eqn. 

the B path are thus 



dip 



K2 

jdkA{k)=-l(^{k2)+wk2) 



(54) 




1 5 10 15 20 25 30 35 40 45 50 
a 

FIG. 4: (Color online) Interpolation between Wannier spec- 
trum (red) and entanglement spectrum (light blue) at the 
edge crossing point with winding number w = 1. N = 50 
unit cells in the a2 direction, with periodic boundary condi- 
tion, is used. Other parameters are the same as those used in 
previous figures: [ti,t2,ts] — [0.3, 0.4, 0.5], m = 1.4, = 0.37r. 
{va} is the spectrum of QRji Q, where a is the level index. 
DiflFerent lines correspond to different values of /3. The red 
dots (/3 — >■ 0) correspond to the limit of Wannier center spec- 
trum, therefore the vertical coordinates are pinned at half 
integers (recall the Berry phase is tt). The light blue dots 
(/3 — > oo) correspond to the entanglement spectrum (rescaled 
and shifted according to eqn. |53[ ), by which the mid gap levels 
at Ta ~ 25.5 translate to an occupancy of / = 0.5, ra = 1 to 
f — and Ta = 50 to / = 1. All curves exhibit particle-hole 
symmetry (with respect to ra = 25.5) due to the coplanarity 
of _B at fcc. Notice the two levels pinned at Va = 25.5 for all /3. 
Had we used open boundary condition in obtaining Q, there 
would be only one such level. 



where (p{k2) = ¥^(^2) — is the non-winding deviation 
of ^p{k2) from a pure winding term. The Wannier state 
corresponding to the I**^ unit cell is a linear combination 
of the Bloch states belonging to the occupied band, 

l*/)=E/fef l*(^2)) , / = l,2,...,A^ (55) 



where |5'(fc2)) is the Bloch state defined in eqn. 42 and 
the coefhcients / ar^iSl 



Jk ^ 



/N 



exp i^—ilk — I f{k)^ 



(56) 



Recall that in the rotated frame, the internal-space Bloch 
states of the occupied band are |'0(fc2)) = ^ (e*'^ )' ^^^^ 
the Wannier wavefunction of the rotated A and B "sub- 
lattices" - which are linear combinations of the original 



A and B sublattices in the same unit cell - are 



1$ 



TiE/if Ifc2> = -;i^^n^)|/), (57) 



I I ~ \/2 ^-"=2 1"-^/ - /S'- 



k2 



(58) 



Here, |/) = e-'-'''^\k2) /VN is the 7*^ reabspace basis 
vector; the unitary profile operator S^ip) is 



^(^)^^|fc2)exp(i^(A:2)) (fc2 



(59) 



and its effect on |/) is to generate a wavepacket centered 
at the /, with its profile determined by the details of 
the fluctuation <f{k). We can then conclude that at the 
coplanar kc points, the two sublattices of the Wannier 
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w= 1 



M+1 M 



FIG. 5: (Color online) Schematics of localization of the two 
sublattices at the edge crossing point with winding number 
w = 1. N = 2M is the number of unit cells. Each empty dot 
represents a unit cell. Red and blue dots represent the A and 
B sublattice wavefunctions of a Wannier state, respectively. 
The A sublattice is w = 1 unit cell ahead of the B sublattice. 
The Wannier state itself is represented as the box around the 
colored dots. The complete set of Wannier states is generated 
by shifting a given Wannier state by integer number of unit 
cells. Since the two sublattices are localized on neighboring 
unit cells, the Wannier center, which is their average spatial 
location, is a half odd-integer. The entanglement eigenstates 
are qualitatively the same as the Wannier states in terms of 
sublattice localization, although quantitatively the wavefunc- 
tions do have projections onto neighboring Wannier states. If 
the bipartition cut passes through the box (representing the 
dominant Wannier state projection of the entanglement eigen- 
state), then each half has only one sublattice manifesting as a 
/ = I occupancy mode. For periodic systems, there are two 
such modes, illustrated in the left and right panels. Clearly, 
only the right panel survives as entanglement f ~ ^ mode for 
open boundary systems. The middle panel is a Wannier state 
whose localization does not cross the bipartition cut. Qual- 
itatively, it represents an entanglement state close to / = 
for the left half system and / = 1 for the right half system, 
cf. the light blue line in Fig. [4] 



states (in the rotated frame) are localized in different 
unit cells with a separation given by the winding num- 
ber w of the path of B. The Wannier centers are spatial 
averages of the locations of the two sublattice wavefunc- 
tions, Xi = I — w/2, thus for odd w, they are half-odd- 
integers. The double degeneracy as evidenced in Fig. |4] 
is due to the periodic boundary condition, because for 
I~w — 1, I — w is N instead of zero, thus its Wannier 
center is at {N + l)/2 which is degenerate with the level 
of / = iV/2-|-l. This is schematically shown in Fig. 5] As 
/? deviates from zero, the eigenstates of GRp Q leak into 
neighboring Wannier states but are still dominated by 
the /3 = component, thus the physical picture stays the 
same. When /3 is infinite, the two degenerate Wannier 
states become the entanglement half occupancy modes 
(before projecting onto half systems) since the entangle- 
ment cut assigns the two sublattice wavepackets into dif- 
ferent half systems, leaving each half system with a "half 
electron" . 

There is an interesting consequence in entanglement 
spectrum due to the w-separation between the sublat- 
tices. For a general \w\ ^ 1, it is obvious from Fig. [s] 
that for each Wannier state, there are w ways to place 
an entanglement cut that will assign the two sublattice 



wave packets into different halves, generating a half oc- 
cupancy mode. Thus in the thermodynamic limit, there 
are 2w modes of / = ^ with a periodic-boundary system, 
and w such modes with an open-boundary system. We 
have confirmed this with hand-crafted (p{k2) relations, 
e.g., f{k2) — w\/2TTk2. This agrees with the general ob- 
servation that the number of entanglement spectral flows 
is given by the Chern index. 



2.3. Armchair Edge 

The armchair edge is more conveniently described by 
a new set of primitive vectors a/: 



Qi = fli 



02 



02 



(60) 



and the associated Bloch phases k[ — k- d[ are 

fc'i = fci + k2 , k'2 = k2. (61) 

The armchair edge is parallel to a[ and the corresponding 
open boundary Hamiltonian can be parameterized by k'^ . 
The bulk Hamiltonian is still eqn. [3| When k[=0, then 
By — and the path of B upon varying ^2 is on the 
zx plane. Thus in the topological phase (origin inside 
the path), one can repeat the same analysis as done for 
the zigzag edge and conclude that at k[ — 0, the Berry 
phase is tt, and the entanglement occupancy is |. We 
have confirmed these numerically. As for the energy edge 
modes, while there appear to be no explicit solution, we 
have confirmed numerically that they still cross at this 
coplanar point. 

From eqn. |3j By /{I + B^) = tan(/c5^/2), i.e., the pro- 
jection of B onto the xy plane is a line passing through 
the point [B^, By) = (—1,0), so the path of B is always 
coplanar for any fixed k'l, but only the one with k'l — 
has the origin on the same plane. Thus unlike the zigzag 
edge case, there exists no other k'^ points whose Berry 
phase can be related solely to the winding of a planar B 
path. 



3. BILAYER HALDANE MODEL 

A bilayer extension can be constructed by Bernal 
stacking two monolayers with a vertical interlayer hop- 
ping i_L, as shown in Fig. [6] t±_ couples any A site 
with the D site above it. The Haldane phases of the 
two layers, for AB and x for CD, are taken to be 
independent. Semenoff masses for the four sites are 
•niA — -^riiB = —mc = Tnjj — m. Both layers have 
the same second neighbor hoppings ti,i = 1,2,3. The 
bulk Hamiltonian is thus 



T = 









(62) 
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FIG. 6: (Color online) Bilayer model. A,_B, C and D label 
sublattice sites, where A and B are on the first layer (solid 
line lattice) and C and D on the second (dashed line lattice). 
Box represents a unit cell. Interlayer hopping is between A 
and the D on top of it. 



where Hab and Hcd are given by eqns. [2] and |3] with 
different parameters Hab = H{(p,m) and TJc^i = 
H{x,—m). Superscript t denotes matrix transposition. 
The two layers individually can have different Chern in- 
dices. Thus as long as the central gap does not close 
when t± is turned on, the total Chern number should be 
given by the sum of those of the individual layers. This 
is numerically verified, see Fig. [7j where we plot open 
boundary energy spectrum and entanglement occupancy 
spectrum. 

An interesting case is when = -X (Fig. [7(d)] ). The 
bulk Hamiltonian is mapped to its complex conjugate 
under inversion, 



IH{k)I = H*{k) , 1 = 



I 1\ 
1 

1 

Vi ) 



(63) 



We shall refer to this as the I* symmetry. If both layers 
are individually topological, their Chern numbers are op- 
posite and there is no gapless edge mode. The entangle- 
ment occupancy spectrum also exhibits no spectral flow, 
but unlike the generic ^ —x case, there are protected 
half occupancy modes occuring at the k\ very close to 
the edge crossing fcc point of the corresponding mono- 
layer (clearly, they must be identical when = 0) . This 
is reminiscent of the Z2 inversion-symmetric topological 
insulators^^, which also shows protected j —\ entangle- 
ment modes without gapless edge modes in the energy 
spectrum. An essential difference is that there, they only 
occur at inversion-symmetric fci = or tt points because 
inversion relates H{k) to H{—k). 

Given the connection between entanglement spectrum 
and Wannier centers, the question naturally arises as to 
how the latter would behave in the presence of T* sym- 
metry. We shall address this question in this section. 



3.1. Non-Abelian Wannier centers and Wilson loop 
phases 

The Haldane bilayer has v = 2 bands filled at half fill- 
ing. For adiabatic evolution of v bands, the concept of 
Berry connection ( "geometrical vector potential" ) gener- 
alizes to a, V X V Berry connection matrix ( "non- Abelian 
gauge field" 



(64) 



where {ip^) is the internal-space eigenstate at k of the 
a*** band, etc. The Wilson loop, also a. x matrix, is 
defined as 



yV = Pexp{i dkA{k) 



(65) 



where P denotes path ordering. The phases of its eigen- 
values play the role of Berry phase. 

The Wannier states and Wannier centers can still 
be defined as the eigenstates and eigenvalues of the 
band-projected position operator eqn. |38[ where Q now 
consists of v bands, Q = X]a=i Sfc ^^^^ 
l*fe) — \k) ® IV'fe) is eigenstate of the a"^ band be- 
low Fermi level. One then finds that as the number of 
lattice sites — > 00, the Wannier centers are given bjff^ 



27r 



(66) 



where / is an integer labeling the corresponding unit cell, 
and 7^ is given by the phase of the w**^ eigenvalue of VV, 



1,2, 



(67) 



The V bands recombine to yield v Wannier functions as- 
sociated with each unit cell. When v = \, this reduces 
to eqn. [39) 

In the finite N case, it is more convenient to replace the 
position operator R — Y®1 with the fc-space translation 
operator eyiY>{iY5k) ® I where 5k = ^tt/N is the fc-space 
step ^^. T he object of concern is instead (cf. derivation of 
eqn.IC22|) 



W = R P, 



P, 



(68) 



where fc„ — ndk with n = l,2,...,N, and Pk^ — 
J2a=i IV'fc )('0fc I is the internal-space projector onto the 
occupied bands at k. Note that this object is basis- 
independent — in particular, it is oblivious to the phase 
choice of the 1-0^) states, which is advantageous for nu- 
merics. Note also that the dimension of Pk and hence of 
W is the number of total bands q whereas that of W is 
v, the number of occupied bands. Ho wever. W is essen- 
tially the nonzero block of W (cf. eqn. C23l, hence they 



have the same (complex) nonzero eigenvalues {pwS^''™} 
where w = 1,2,...,;/. The amplitudes will deviate 
from unity due to discretization. The phases 7^, will be 
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FIG. 7: Open boundary energy spectrum (upper panel) vs. entanglement occupancy spectrum (lower panel) of the bilayer 
model. Gray line indicates Fermi level. Parameters used are [ti,t2,t3] = [0.3,0.4,0.5], m = 1, and t_i_ = 0.5. The Haldane 
phases </!> and x for each case are given in the subfigure title, (a) Both layers have Chern number C = 1, giving a total "^0 — 2. 
The entanglement occupancy spectrum has two flowing branches, each taking about one fei cycle (fci — >■ fci + 27r) to flow from 
/ ~ 1 to / ~ 0. (b) AB layer has C = 1 and CD layer has C = 0, hense "^C — 1. Note the avoided crossing near / ~ 0.1: 
If the two Haldane layers are not coupled, the crossing will not be avoided and one ends up having one flowing branch (which 
takes one ki cycle to flow from / ~ 1 to / ~ 0, see (a)) and another trivial branch. The avoided crossing effectively merges the 
two branches into a single branch which now takes about two fci cycles to flow from / ~ 1 to / ~ 0. (c) and (d): AB and CD 
individually have opposite Chern indices, (c) Generic case where (p and x a-re not related, (d) I*-symmetric case with (p = —x- 
Insets: details around half occupancy / = |. In (d), the crossing at / = | is not lifted by the t± hopping. Crossings near 
/ = 1 and are avoided in both (c) and (d). Note that in all cases the number of gapless edge modes at a given boundary is 
given by C, so is the number of entanglement spectral flows, (d) has no flow due to avoided crossing near / = and 1. 



interpreted as the Wannier center offsets as in eqn. 66 
See Appendix [C| for details. 
W in the form of eqn. 
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as a product of occupied 
band projectors over the period of fc - is also known 
as a monodromy and has been used in analyzing e.g. 
inversion-symmetric Z2 topological insulators'^. One 
question is whether or not it matters to break up the 
Wilson loop at some k point other than fc = 0. The an- 
swer is no. This is because the eigenvalues of a product 
of (possibly non-commuting) projectors are invariant un- 
der cyclic permution of these projectors. We prove it in 
Appendix [P} 



3.2. 



Wilson loop phase of the X*-symmetric bilayer 
Haldane model 



3.2.1. X* -symmetric 



r eigenstates 



Upon swapping the C and D rows and columns, the 
Hamiltonian eqn. 62 becomes (suppressing fci depen- 



dence) 



M T 

T* M* 



(69) 



M = B -a 



T = T* 



t_L 




(70) 



where uj and B are given by eqn. [s] The form of eqn. 
means its eigenstates are of the form 



69 



u*' 



(71) 



where \u) = {ul) is a two-component column vector. 
The eigenvalue equation is thus 



M\u) 



-T\u*) = e\u) 



(72) 



Since interlayer hopping is restricted to A and D sites, 
T\u*) only depend on wi, thus one may always adopt 
a phase choice of \u) such that Ui is either real, in 
which case T\u*) = T|m), or imaginary, in which case 
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(a) (f> = O.Svr, X = -OAtt 



(b) <f> = O.Svr, X = -O.Stt 



FIG. 8: (Color online) Non-abelian Wilson loop phases. Haldane phases are given in subfigure titles. Other parameters used 
are the same as in Fig. [7] Insets: details around 7 = vr and 2tt. (a) has avoided crossings near both j — tt and 2-k, but (b) has 
no avoided crossings. The red dots in (b) can be computed by eqn. 80 The blue dots are given by their negatives (shifted by 
2n). 



T\u*) = —T\u). The eigenvalues of H consequently com- 
prise those of 2 X 2 defined as 



=uj] 



M ±T = uj±I + 
UJ± = uj ± ltj_ , 



B± ■ a 



(73) 
(74) 
(75) 



where (|B±|,z?± , are the spherical polar coordinates 
of B±. The role of t± is two-fold: by modifying u), it 
splits the two monolayer copies; by modifying Bz, it also 
changes the level splitting of each monolayer. Thus as 
t± is increased from adiabatically, it is possible to re- 
arrange the order of the monolayer bands. Here we shall 
focus on the case where the central gap is never closed 
during — >■ t_L, so that the occupied bands \ip^{k)) at 
half filling are generated by the isospin-down states of 
, viz.. 



72 



(76) 



where Dz{f) 



diag(l,e*^,l,e-^'^) and \'d^ 



cos 



are purely real. Note that have the 



prescribed form 
real for IV" )■ 



with ui imaginary for and 



3.2.2. Wilson loop phases 



Fig. |8] plots the Wilson loop phases of the Haldane bi- 
layer where the two layers have opposite Chern indices. 
In the generic </> ^ — x case (Fig. 8(a) I, both 7^ lev- 
els exhibit no fiow over the period of ki due to avoided 
crossings, which is similar to the corresponding entan- 
glement spectrum (Fig. 7(c)). In the I*-symmetric case 
where 4> — — X: none of the crossings at 7 = and tt 
are avoided, thus both 7^ levels flow in opposite direc- 
tions .^This is different from the entanglement spectrum 
(Fig. 7(d) ) where the crossing at / = ^ is preserved, but 
the crossings near f — I and are avoided. Unlike the 
monolayer case, the f — \ crossing does not coincide in 



ki with the 7 = tt crossing, although they are very close. 
This is because one no longer has the analog of a coplanar 
ki point. We will come back to this point later. 

Since the set {7^,} is periodic as fci —> /ci + 27r, their 
winding numbers \jw{2-k) — 7.^,(0)] /(27r) must be inte- 
gers. The question is how to extract them. From eqn. |76[ 
one finds that the Berry connection matrix is purely off- 
diagonal. 



(V'+|\|V'+) = (V'-|\I^-) 



0, (77) 

cos(ii?+)cos(ii?-) , (78) 
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hence 



W(/ci) =e 

V(27r) 



i'y(ki) ay 



dip COS ^ COS ^ 



(79) 
(80) 



v(o) 



Note that when tj_ = Q, d'^ = d~ and eqn. 80 reduces 
to e qn. [45l The Wilson loop phases are simply ±7. In 
Fig. 8(b) red curve corresponds to +7 and blue to —7. 



While there are no inversion or time-reversal symme- 
tries protecting the winding of this off-diagonal 7 (by 
fixing 7 = and tt at T or I symmetric fci points), it 
is nontheless robust with respect to as long as the 
central gap remains open, and is thus given by the cor- 
responding monolayer Chern number (tj^ = case) . To 
see this, we use the eigenstates of VV (cf. eqns. 
C14| to construct a set of states. 



C13 



and 



9l IV^+> 



9± \^-) -- 



1 

72 



{|^+)±z|V'-)} (8 



1) 



(82) 



where g± = {g'^^g^)^ are the eigenstates of VV, in this 
case, the eigenstates of ay. The Berry connection matrix 
A is diagonal in the \ii±) basis, i.e., 7 is the Berry phase 
(over of the I77+) "band". It is easy to verify that the 
two layers, AB and CD, each contributes exactly one 
half to 7. One may thus construct a ficticious two-band 
model whose occupied band is given by the projection of 
say 1 77+) onto the AB layer, 



1 

72 



1 



and the unoccupied band as 
'l 



Here, 



1+) 



/ c 



1 

\/2 



(83) 



(84) 



■cos(z9/2)\ 
Vsin(T?/2)j 
of By construction, { 
tonian is 



are the iso-spin up counterpart 
+ 1— ) = 0. The ficticious Hamil- 



H = 



es{k)\k){k\®\s{k)){s{k)\ (85) 

k.s=± 

where e+ can be chosen as (say) the lower energy of the 
two unoccupied bilayer bands, and e_ as the higher en- 
ergy of the two occupied bilayer bands. Then 7(^1) is 
the Berry phase of H, and its winding number is the 
same as the monolayer Haldane model as long as t± does 
not collapse the bilayer central gap. The reason why the 
point of 7 = TT is very close to but no longer coincide 
with where the / = ^ entanglement mode occurs, which 
was the case for monolayers, is that if one maps the state 



of eqn. 83 back to a vector B, with its azimuth as Lp and 
its polar angle being some sort of average of and , 
then around the monolayer kc, the path of B comes close 
to coplanarity (yet never exactly so) . Since the "temper- 
ature" j3 is different for the entanglement and Wannier 
spectra (cf. eqn. 51), they need different distortions in 
the B path to balance f3 in order to reach their special 
values. 



3.3. Generic X*-preserving interlayer hopping 

The reason why 7 can be extracted from the bilayer 
Haldane model is because the Berry connections A{k) 
are proportional to ay, and hence are mutually commut- 
ing at different k. This in turn is because the interlayer 
hopping is only between A and D sites, so that H de- 
composes into sectors. Thus while there is a ma- 
trix structure, the situation is nontheless Abelian. Still, 
the I*-symmetric bilayer provides an interesting exam- 
ple where the topological information is not obvious from 
the open-boundary energy spectrum. The role of the I* 
symmetry is in prescribing the eigenstates to a particu- 
lar form ( j^^l^ ) (after swapping C and D labeling). One 
question to ask is if the crossings as seen in the Wilson 
loop phases of the I* -symmetric bilayer is a consequence 
solely of the I* symmetry, or if it also depends on the 
interlayer hopping being restricted to ij^ only. The most 
general interlayer hopping that preserves X* is 



T = 



Ai A2 

^2 A3 



A, € 



(86) 



instead of the one used in eqn. [70] In the bilayer stacking, 
the B sites are surrounded by three "nearest-neighbor" 
C sites on the other layer, C by A, and D by B, thus 
one example of such T matrix is by treating the strength 
of all these interlayer hoppings as A in addition to the 
vertical t± between A and D, 



A2 - A(l + e*^-i +6''=^) , 
A3 = \{e'^^ _^g»fc2 _^g»(fei+fe2)) 



(87) 



Although not important for our purpose here, these 
parameters can be connected to the SWMC model of 
graphen^^. 

I* symmetry of the Hamiltonian is inherited by the 
internal-space projector of occupied bands, and hence 
their product W, 



IP{k)I = P*{k) , IW(fci)I = W*(A:i) 



So if \gyj) is an eigenstate of W with eigenvalue 



' , then T \gw)* is also an eigenstate with eigenvalue 
Thus the eigenvalues of W come in complex- 



Pu 
Pu 

conjugate pairs. For the bilayer with a general interlayer 



15 




(c) <t> = 0.3n, X = -OAn, N2 = 10 (d) tj) = O.Stt, x = -O.Stt, N2 = 10 

FIG. 9: Wilson loop phases of bilayer Haldane model with generic I*-preserving interlayer hopping, and the eflFect of k2 
undersampling. Parameters used are t± = 0.5, A — 0.4, other parameters are the same as in Fig. [7] Insets: details around 
7 = TT and 2-K. (a): 7^ —x, i.e., I* -non-symmetric, (b): (j) = — x, i.e., I*-symmetric. (c): same as (a) but only use N2 = 10 
projectors in eqn. |68[ (d): same as (b) but with = 10. Notice that qualitative features do not depend on the number of ^2 
discretization. In particular, k2 undersampling does not lift the crossing in the I*-symmetric case. 



hopping eqn. |86[ while there does not seem to be a way 
of tracking one of the two W eigenstates (with nonzero 
Pu,) analytically, we do find numerically that the {7^1} 
pair flow in opposite directions and cross at 7 = and 
TT just like the case with t±-on\y interlayer hopping, and 
that their winding numbers are given by the underlying 
monolayer as long as the bilayer central gap does not 
collapse when turning on tj_ and A. Fig. |9(b)] plots a I*- 
symmetric case with T given by t± — 0.5 and A = 0.4. 
Compare with Fig. 9(a) where <j> ^ — x and the flows are 



broken near 7 = tt and 2tt. 

An interesting observation made possible by express- 
ing the Wilson loop as a product of projectors is how the 
topological and non-topological cases behave with an un- 
dersampling of , that is, when the number of projectors 
used in eqn. |68] is reduced. For example, if one were to 
reduce N projectors to say N/3, it is equivalent to replac- 
ing every two out of three projectors in the monodromy 
by a unity operator, thereby allowing wavefunctions to 
leak into unoccupied states when "propagating" from ki 



16 



to ki-3. Of course W ceases to be unitary due to the 
leakage, which is just the reason why < 1 in the fi- 
nite case. A somewhat unexpected behavior is that 
the Wilson loop phases in the I* -symmetric case now 
become degenerate for an extended region of fci values at 
J — TT and 0, as opposed to crossing at discrete points in 
the N ^ oo limit. The eigenvalues there are simply real 
numbers, which are their own conjugates and thus not 
ruled out by the I* symmetry. Note however that the two 
real solutions are not required to be the same — indeed 
they are different for finite N and only approaches ±1 as 
N ^ oo per unitarity of VV. No such phase degeneracy 
is observed for the generic I*-asymmetric cases, where 
the phase flows still exhibit avoided crossing. The topo- 
logical signature is in this sense more prominent away 
from the thermodynamic limit, unlike e.g. the crossing 
of the monolayer energy edge states which in general is 
a thermodynamic limit result. 



3.4. Relation with a Z2 inversion-symmetric TI 

Hughes et. al. studied a Z2 inversion-symmetric TI 
model (HPB model)i2i, 

H{k) =smk^Ti+ sin kyT2 + M{k)Ta + B^Tb , (89) 
M{k) = 2 — m — cos k^ — cos ky , (90) 

Tl=(J^®T^ , T2=l®Ty , (91) 

ro=I(^r, , Tb^o,®{\) . (92) 

Here, both ai and are the Pauli matrices, cr, act on 
the spin indices and Ti on the orbital indices. The HPB 
model is built from the BHZ quantum spin Hall modeP^l 
by adding the term which breaks time- reversal but 
preserves inversion: the inversion operation is defined as 
T' = Fo such that VHik)!^ = H{-k). The authors 
showed that it has a Z2 topological index protected by 
the 7^-symmetry. 

The HPB model is also I*-symmetric: I — ® 
such that IH{k)I — H*{k). In fact this is the mirror 
symmetry as noted by the same authors in Ref. 1121 The 
T-breaking term plays the role of t± of the bilayer 
model, therefore the model can be cast into the form of 
eqn. 69 by rotating {t^ 



,Ty) to {Ty,—Tx). One can then 
solve it analytically and analyze its Wannier centers in 
the same way as the bilayer model. Its Z2 index is equiv- 
alent to the winding number of 7, and is inherited from 
the "monolayer" (individual spin species) as long as the 
central gap does not collapse. 

It turns out for this particular model, one can choose 
to break either V or I* (but not both) and still retain 
the nontrivial topology. Here we only consider breaking 
the V, for as long as it is preserved, the proof of Ref. [T2| 
applies. To break V, we add in an I* -symmetric "spin- 
fiipping" hopping. The most general form is to replace 




FIG. 10: I*-symmetric HPB model with broken inversion. 
Top panel: periodic boundary entanglement quasi-energy. 
Bottom panel: Wilson loop phases. Parameters: m = 1.1, 
(Ai,A2,A3) = (0.7,0.5,0). Since the I* symmetry is pre- 
served, the model still retains its nontrivial topology: in the 
entanglement spectrum, although there is no spectral flow, 
the £ = mode, which now is slightly shifted away from 
kx = 0, is still robust with no avoided crossing; for the Wilson 
loop phases, the two branches still wind in opposite directions 
and intersect at 7 = and tt, but the crossings are no longer 
pinned at fc^, = or tt. Insets: details near y — tt and 2tt. 



the Bx Tb term in eqn. |89| by 
T' 

Q I ' 



Al A2 

-A2 A3 



(93) 



are complex numbers which could have 
The extra "— " sign in T as compared 



where again Ai 
k dependence, 
with the bilayer case (eqn. 86) will drop after the afor- 
mentioned (r^jTy) — ^ (tj,, —t^) rotation. One can easily 
verify that the condition for inversion symmetry to also 
hold is (Al, A2, As)^^ = (Ai, -A2, Aa)^. 
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Fig. 10 plots the entanglement spectrum and Wilson 
loop phases of an X* -symmetric HPB model with broken 
inversion symmetry, m — 1.1, (Ai,A2,A3) — (0.7,0.5,0). 
The topological signatures still persist: the entanglement 
quasi-energy spectrum exhibit robust zero modes, and 
the Wilson loop phases flow in opposite directions and 
cross at 7 = and tt. Since the A2 term explicitly breaks 
V 1 the kx points where they occur are no longer pinned 
at symmetry-invariant points. 



4. SUMMARY 

We have studied the monolayer and bilayer Haldane 
models and identified several topological signatures from 
the real space perspective. The monolayer zigzag edge 
modes can be analytically solved using the Ansatz that 
the wavefunctions of its A and B sublattices are pro- 
portional. This particular form poses restrictions on the 
boundary condition which can no longer be prescribed 
(as e.g. an open boundary) but must be solved self- 
consistenly, and the edge state is recognized as an open 
boundary one when the "tunnelling strength" p of the 
two boundaries approaches zero in the thermodynamic 
limit. Using the edge solution, the transverse momen- 
tum ki — kc at which the two edge modes cross can be 
identified as the coplanar point of the B vector that gen- 
erates the bulk Hamiltonian H = ui + B ■ a — that is, at 
this point, varying the longitudinal momentum ^2 from 
to 27T will drive i? in a closed path on a plane that 
passes through the origin. The problem of mapping the 
2-D Brillouin zone (fci, k2) to the sphere B{ki, ^2) is thus 
reduced to a 1-D mapping from fc2 to the ring B{kc, k2), 
and the bulk Chern index reduces to the winding number 
of the B ring around the origin. 

Interestingly, at fee, both the entanglement spectrum 
and the Wannier centers exhibit crossings similar to the 
edge modes. If one treats the entanglement spectrum as 
the bipartition coarse graining of the Wannier centers, 
then there is a continuous parameterization /3 such that 
the eigenvalues {ra} of a band-projected real space op- 
erator QRp Q correspond to the entanglement spectrum 
at one limit /3 — > 00, and to the Wannier centers at the 
other limit /3 — )■ 0. The crossing at fee is universal in 
that for arbitrary /?, there always exist levels fixed at 
Ta = {N+1) /2, where N is the number of unit cells in the 
longitudinal direction 0,2 . The special value translates 
to half occupancy / = ^ in the entanglement case, and to 
a state localized right in the middle of the iV-cell chain in 
the Wannier case. This universal crossing can be traced 
back to the B coplanarity because the two sublattices 
(after an internal rotation) are localized in different unit 
cells with the separation given by the winding number 
w of the B ring, thus their average center is a half-odd- 
integer if w is odd, giving rise to the special Wannier 
center, and an entanglement cut placed in between them 
will assign the two sublattices to different halves of the 



bipartite, yielding an / = ^ entanglement occupancy in 
each half. 

The Haldane bilayer is constructed by Bernal-stacking 
two monolayers and allow vertical interlayer hopping i j_ . 
Without exception, the Chern index of bilayer at half 
filling is the sum of those of the individual monolay- 
ers, which also manifests as the number of entanglement 
spectral fiows. A special case is when the two mono- 
layers individually have nonvanishing Chern indices but 
are opposite to each other: the nontrivial topology sur- 
vives in a sense if their parameters are exactly opposite 
such that the bilayer bulk Hamiltonian is mapped to its 
complex conjugate by inversion, XH{k)I = H*{k) {I* 
symmetry): while there is no gapless edge modes and no 
entanglement spectral flow, the entanglement spectrum 
does exhibit protected f = \ modes. It becomes more 
prominent in the Wannier center (Wilson loop phases) 
spectrum, where two branches start to flow in opposite 
directions with the magnitude of their winding number 
given by the monolayer Chern index. The flow is ro- 
bust as long as the central gap of the bilayer does not 
collapse when t±^ is varied. The role of X* symmetry is 
further confirmed numerically by adding in I* -preserving 
generic interlayer hopping. We also found that unlike the 
monolayer open boundary edge spectrum whose gapless 
crossing is a thermodynamic limit result, the crossing of 
the opposing Wilson loop phase flows are more promi- 
nent as the number of unit cells N in the 0,2 direction is 
reduced: they are now degenerate for an extended region 
of ki instead of crossing at discrete points. This does 
not happen for the topologically trivial case where T* is 
broken. 

The flow-less entanglement spectrum with protected 
midgap modes is reminiscent of the Z2 inversion- 
symmetric topological insulators, with the essential dif- 
ference that the midgap modes there are pinned at 
inversion-invariant k points. We looked at one such 
model studied in Ref. T^, This particular model has 
both inversion {V) and mirror symmetry. The latter co- 
incides mathematically with the X* symmetry, rendering 
the model solvable in the same way as the t_L-only bi- 
layer. The nontrivial topology as captured by the Z2 
index survives as long as either of V and I* is preserved. 
The former case falls within the analysis of Ref. [T^ We 
illustrated the latter by adding in I* -preserving pertur- 
bations that break V. Such a perturbations shift the 
f = \ modes and the tt Wilson loop phases away from the 
■P-invariant k points but does not generate any avoided 
crossing. 

Both the BHZ quantum spin Hall model and the HPB 
inversion-symmetric TI model inherit their Z2 index from 
the Chern index of the underlying single spin species, 
and from the point of view of the Wilson loop phases, 
are topologically equivalent to the situations where the 
coupling between the two spin species are turned off, al- 
though for the HPB model this qualitatively changes the 
edge spectrum behavior. Mathematically the question 
becomes how the interspin coupling (or interlayer cou- 
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pling in the bilayer case) influences the flow pattern of which can also be checked explicitly using eqn. 23 Notice 
the Wilson loop phases, and what kind of coupling does 
not degenerate the pattern to the trivial case {e.g., that 
of a unity operator). The bilayer Haldane model shows 
that there are other types of coupling which break both 
time-reversal and inversion and yet still result in a non- 
trivial topology. 



that this applies to the singular case, eqn. 26 as well. 

That we can get eqn. |A7| is actually fortunate, for oth- 
erwise there will be no twisted boundary consistent with 
the Ansatz. We will come back to this point later when 
deriving the eigenstates. 
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Appendix A: Details of Monolayer Zigzag Edge 
States 



1. Edge solution 

From [16} we note that as fci — — fci, 
wi(-fci) = ^2(^1) - Pi(-fci) 



■Pi{ki 



(Al) 



Eqn. 21 gives two equations (the ratio r itself being yet 



undetermined) with two unknowns, A and e, for each ki 
value. The second equality of eqn. [21] does not involve e, 
and can be used to solve for A, 



P2 ^2^^ + (^^2^1 - '^i^^2 + P2) ^ + P2 «i = 0, 



(A2) 



This yields eqns. 22 and 23 Note that sending A O A^^ 
and vi o V2 keeps eqn. A2 invariant, thus together with 
eqn. |A1[ one gets 

d{k) = d{^k), (A3) 
\Z^{-k), r+{k) ^rZ\-k). (A4) 



A+(fc) 



To solve for e, we use a simple fact about ratios: if 
r = xi/yi = X2/y2, then r is preserved by arbitrary 
linear combination of numerators and denominators, r — 
[axi -t- bx2)/{ayi + 62/2), except the unfortunate choice 
that makes axi + bx2 — 0. Applying to eqn. 21 we get 



P2/11 -P2£ - 2piRez)i 

r = — , — (AS) 

P2h2 - P2S - 2piRev2 

where Re indicates real part. This gives the two solutions 



to eqn. |A3 



of e, eqn. 24 in terms of r± as solved in eqn. 23 Similar 



-(m, fci) = £_(— TO, —ki). 



(A6) 



The bulk spectrum has the same inversion property 
which is easily verified from eqn. [3[ 

In the text (eqn. 27) we mentioned that r± is real if 



the discriminant A is non-negative. The reality of r± 
(eqn. 271 in turn has the following implication: Using 



r — r* and eqn. 21 



X*v* 



X*V2 +P2 



r = -|Ap, (A7) 



To derive the edge crossing point, we note that ac- 
cording to eqn. [23] the ratio r depends on the branch: in 
general r_|_ ^ r_. But eqn. X5 implies that at the edge 
crossing point (s), ~ r^. It would seem the latter is 
satisfied only when the discriminant A = 0, but one can 
easily verify from Fig. [2] that these do not correspond to 
the edge crossing points. In fact, there is a range of pa- 
rameters in the topological phase where A is always posi- 
tive, e.g., ti = t2 = — 0.3, TO — (p/i: — 0.5. Recall that 
in deriving eqn. |A5[ we used linear combinations of de- 
nominators and numerators of the three ratios in eqn. |21[ 
the validity of which requires these combinations to be 
non-singular (cf. discussion leading to eqn A5 ) . Thus the 
only way for eqns. [23[ and [X5[ to be consistent — the for- 
mer implying r_|_ ^ r_ , the latter implying otherwise — is 
for the linear combinations of both the denominators and 
the numerators to be singular, viz.. 



P2 



{hi - e) - 2piRei;i = P2(^2 - e) 



2piRew2 = 0. 

(A8) 



This yields the edge-crossing condition eqn. 



3. Eigenstates 

Eqns. [19] and [20] can be cast into the same Schrodinger 
equation. 



(A9) 



where a, h and c can either be the set of numerators or 
denominators in eqn. |21[ e.g., a = vi, b = hi + Xpi — e, 
c = v\ + Ap2- Note that these are already known once 
the edge solutions are obtained. Eqn. |A9|is equivalent to 



V'„+l - Xil/'n = X2 (V'n - aJl-^n-l) , (AlO) 

(All) 



b c 

Xi+ X2^ , XiX2^ - 

a a 



i.e., Xi and X2 are solutions to 

ax^ -|- 62; -|- c = 



(A12) 
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Denote (/>„ = 
X2~^4>i, and 



V-'n — a^iV'n-ii then by eqn AlO 



2^1 </>n-l + <pn-2 ' 
„n— 2 J I _n— 1 



(1) If xi X2, the geometric series can be summed, 

V'n = — — 01 + Xi ■00 = /n '01 - a;i a;2 /„-i'0o , 



Xi — X2 



/n — 



(AM) 
(A15) 



(2) If xi = X2 = X, then 

V'n = nx"- Vi - (?^ - l)a;"Vo , (AI6) 

which is the same as applying L'Hospital's rule on the 
previous case. 

(3) If ac = 0, one of the x, say X2, is 0, then 



V'n = X'l Vl ■ 



(A17) 



Thus in all cases we may proceed with eqn. |A14[ This 
yields 

i'N+i = /w+101 - X^X^fjy^l^O , (A18) 
We now come to the issue of boundary conditions. 



Eqn. 17 implies 



N 



V A0^+i / VAV-i 



From eqn. |9j 



1 



V2 



W1U2 V~P2 Vl 



(A20) 
(A21) 

(A22) 



which exists if W1W2 7^ 0. The case V1V2 — can be ana- 
lyzed by Taylor expanding in the vanishing Vi similar to 
the discussion leading to eqn. 26 The can be dropped 
from eqn. |A20"1 For eqn. |A2l"j one gets from eqn. [2T] that 



Vl r 
7 U 



(A23) 



thus the two boundary conditions eqns. A20 and A21 
become 

"00 = (]^ , 0Ar+l = pAU . 

(A24) 



(A13) T^i^^s implies both ( \ ) and ( x ) are eigenstates of U, 



(A25) 



Thus, iiU ^ e*^I, then ( \ ) and ( ^ ) are either equivalent 
(if 61 =02), indicating 1 = r, or orthogonal (if 9i 7^ 62 



indicating r — — |Ap. The latter is nothing but eqn. A7 
The boundary conditions then reduce to 



■00 = pe ''^'^JN , -0^+1 = 



(A26) 



Substituting these in eqn. A14| yields the following con- 
dition. 



xi'J:2fN-l ^ije 
/iv + l 



Jn+i 



Tn + 1 



1 = (A27) 



= B 



where one needs to tune 9i and 02 such that at least one 
solution of p is real. A sufficient condition is for both 
A and B to be real. Here we make A > so that the 
solutions of p are always real. 



arg 



7V+1 



XlX2fN--i 



then 

B ^ Bie-^^i 
thus for B to be real, 



A = 



B2e'^e'^' 



XlX2fN-l 
fN+1 



(A28) 
(A29) 



= B -B* 



B^e 



-i5\ -iBi 



h = arg (Bi - B*2e- 



(A30) 
(A31) 



where the freedom { } can be used to switch the sign of 
B (which has been made real). Then we have 



\B\ ^\B[^+4A 
2A 2A 



(A32) 



Since \B\ and A are both > 0, the p with smaller mag- 
nitude (which has a better chance of — >■ to represent 
an open boundary) is always p+ . Notice that each of p± 
still depends on which branch of the edge solution we 
have picked in calculating xi and X2- 
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In particular, one can show that p = 1 and \xi\ = 1 oc- 
cur simuhaneously, the former means periodic boundary 
condition, while the latter means the Ansatz solution is 
a bulk solution, so this makes sense. To see this, assume 



p = 1 and X2 = e**, then eqn. A27 becomes 



1] 4 



N+l , i4, 



gi0gie2 _ gi(Af+l)0 



(A33) 



One then gets e'^^ = e'^^ — e*^"^ by requiring all square 
brackets to be zero, i.e., eqn. A27| is indeed consistent. 



4. Graphene and boron nitride 



For graphene, K'-^'> — K^'^'^ = 0, so eqns. 19 and 
become 
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R = 



/Pi 

P2 Pi 

P2 Pi 

V 



P2 Pi 



Pi = -2 COS I , P2 = -1 



(A34) 
(A35) 
(A36) 



where z controls the boundary condition. 

If A 7^ and 1/A 7^ 0, i.e., both A and B sites have 
charge density, then Rip a = eV'^/A gives 



P2 i'n + Pi -011+1 = -^0n+l , 

i}o = z'4!M, n = 0, 1,2, . . . ,7V 



(A37) 
(A38) 



thus 



0, 



n+l 



P2 



J -Pi 



n+l 



P2 



-N 



whereas R^ipA = X^ipA will give 



Pi Ipn + P2 Ai+1 = >^eipri , 

ipN+i = zipN, n = l,2,...,N,N + 1 



yielding 



Ae - Pi Y , f Xe- pi 

Vi , z = 



i^n+l 

P2 J \ P2 

Equating the expressions for z gives 



N 



A = ±l 



£j-pi 
P2 



N 



(A39) 



(A40) 
(A41) 



(A42) 



(A43) 



where e = Ae = e/A. Equating the two recursions then 
yields 



e = Pi ± P2 = -2 cos I =F 1, 



(±1)'^ = ±1. (A44) 



Since \z\ = 1, these states correspond to bulk states with 
periodic or antiperiodic boundary conditions. 
If A = 0, then 



£ = , Vn+l 



Pi 
P2 



and 



N 



= -(2cos|)^ 



|z|«i fce(f,f) 
U|>1 fce( 



2ir 27r ^ 
3 ' 3 > 



(A45) 

(A46) 

(A47) 
(A48) 



thus it is an edge solution localized on the A sites at 
the y = 1 edge if |fc| > 27r/3 (i.e., connecting two Dirac 
points), with zero energy. 
Similarly, if 1/A = 0, then 



R 

£ = 







= £ 



0n+l = 





_P2 
Pi 



and 



N 



-(2cos|)^ 



|z|«i fce(f,f) 
|z|»i fce(-f,f) 



(A49) 
(A50) 

(A51) 
(A52) 



thus it is an edge solution localized on the B sites at the 
y = N edge if |fc| > 27r/3 with zero energy. 

For boron nitride (BN), the A sublattice is nitrogen 
and B is boron. If A 7^ 0, 1/A 7^ 0, then XRtpA — — 
m)'ipA yields 



P2 



+ 1 



where £_ = — 



£- -Pl 



N 



(A53) 



Pl \ P2 

Similarly R^tpA = X{e + m)^A gives 



n+l 



£+ -Pi 



Pl 



(A54) 



P2 \ P2 

where £+ = A(£ + m). Equating expressions for z gives 



£- -Pi 

P2 



g+ -Pi 

P2 



(A55) 
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while equating the two ■0 recursions gives 



P2 



P2 

e+ ~pi 



thus 



Pi 



P2 



£+ -Pi 

P2 



= ±1 



= ±v^ 



{Pl±P2Y 



(A56) 

(A57) 
(A58) 



When \z\ = 1, these are bulk states. 

If A = 0, then the edge solution is the same as the 
corresponding graphene case with s = m and the edge 
state is purely on A sites. If 1/A = 0, then e = —m and 
the edge state is purely on B sites. 



Appendix B: Proof of entanglement half occupancy 
modes for Zigzag edge with t2 — 

At the edge crossing point, the path of B is coplanar. 
Rotating to that plane, the internal-space part of the 
occupied bulk states are given by eqn. 43 with i3 — 7r/2, 



The occupied band projector is thus 

g = J2\k){k\®Mk)){m\ 

k 

= i^|/c)(/c|®(I-a^(fe)) 



(Bl) 

(B2) 
(B3) 



where cr^ is the spin operator polarized along the direc- 
tion in xy plane with azimuth ip. 



e-'^P 




The restricted correlation matrix is thus 
G= I^Ffc® (I-(T^(fe)) 



k 

1 _ 

2 '2 



where 



Pfe = R\k){k\R^ 



(B4) 

(B5) 
(B6) 

(B7) 



is the Bloch projector |A:)(fc| restricted to half space, and 
i? is an oblong matrix to project out the first half. 



/I 



R 







0/ 



(B8) 



We have used the fact that = I is the unity in 

the half system. Instead of the entanglement spectrum, 
which is the eigenvalues of G, we focus on the spectrum 
of 



(B9) 



An entanglement half occupancy mode / = 5 corre- 
sponds to a zero mode of G. The order of direct product 
does not influence the spectrum, thus G may be written 
in the following off-diagonal form, 



G 



M 



M = RMR^ 



-iipik) I 



(BIO) 
(Bll) 



Assume the full system has 2A^ unit cells, then the 
allowed k are 



mr 



n= 1,2, 



,27V 



(B12) 



We shall cut the system in equal halves. Let R be the 
complement projector to R, 



1 



R 







1 



(B13) 



V V 

Then for the real space Bloch waves |fc„). 



P|fc„) = 



iNkr 



R\kn) = (-l)"i?|fc„) 



Then one may write as 
M - 



M Q 
Q M 



(B14) 



(B15) 



where Q = RMR^ . Then M and Q can be separated 
into even and odd parts. 



M^Me + Mo , Q^Me-Mo 



where 



Me= J2 -^1^" 



{kn\R^ 



Mo^Y. R\kn)e-"^^''-\k,,\R 

n odd 

It is also easy to verify that 
f2ME 



(B16) 

(B17) 
(B18) 



U'^MU 



thus 



2Mc 



det (M) = det (AMeMo) 



(B19) 
(B20) 
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In fact, the two sets {i?|fc„)} with even and odd n are 
both complete sets in the half space, but normalized to 
l/\/2 due to the projection, thus 2Me and 2Mo are both 
unitary. 



MeMI = MqMI 

and 

det {2Me) = e"'*^ 
where 



det {2Mc 



(B21) 



(B22) 



(B23) 



odd 



Now, when t2 — t^, both vi and V2 are real (eqn. 16), 
so at ki — fcr,, c hanging k2 — > — fc2 only changes the sign 
of By (eqn. 34). Thus in the rotated B plane, (p{—k) = 
—ipik). In particular, 

^(tt) = or tt (B24) 
and (/3(0) is defined as 0. Me and Mo are both real, and 
det [AMeMo) = e-*(*^+*o) = g-^'^^'') = ±1 (B25) 
because all other ip{k) are cancelled by (p{—k). Then 

det (MeMe + MqM'^) 



det Af = det (M^ + Mq) = 



det (MoM(^ + MoMl) 



del Me 



det Mo 
det Mb' 



det Mb 
deti\/'^ = 



(B26) 
(B27) 



-detM"^ if^(7r) = 
-detM'^ if .^(tt) = vr 

(B28) 



where in the first line, we multiplied and divided by 
detM^, and in obtaining the second line, we rep laced 



MeMe with MqMq which follows from eqn. 
det G = (det Mf=0 if ^(tt) = vr 



B21 



Thus 
(B29) 



i.e., there exist entanglement half occupancy modes. 



Appendix C: Non-Abelian Wannier centers and 
Wilson loops 

Consider, for simplicity, a periodic 1-D system. 
Higher-dimensional systems can be analyzed by parame- 
terizing the system with momenta (/C|| , kj_) and consider- 
ing the effectively 1-D system at fixed fcj^ . Assume there 
are q bands and N unit cells. The full Hamiltonian is a 
qN X qN matrix H and its Fourier transform is a q x q 



matrix H(k). The eigenstates of % are giV-component 
Bloch states \^%), 



\nj = \kn) ® 

27m 



N 



a = l,2,... 
n = l,2,...,iV 



(CI) 
(C2) 



The g-component iV'fe ) is the a band eigenstate of 
H{k), and the A^-component |fc„) is the Bloch phase, 
{x\kn) — exp{iknx) / \/N . The Wannier states can be 
defined as eigenstates of GR^^, where 



N 



a— 1 n— 1 

is the projector onto the i> occupied bands, and 
R = exp{iXSk) (g) I . 



(C3) 



(C4) 



X is the position operator in real-space (i.e., measuring 
unit cell coordinates), and Sk — 2tt/N is the step in the 
discrete {fc„}. Note that exp{iXSk) is the momentum 
space translation operator. 



exp{iXSk)\k„ 



Sk) 



(C5) 



this ensures the resulting Wannier states to have proper 
translational symmetry in real space. The eigenstates of 
GRG belong to the occupied bands and thus have the 
decomposition. 



N 



a—1 71—1 



(C6) 



To solve for it is useful to introduce the v x i/ 

overlap matrix Umn between the internal-space bases at 
different km and A;„ points. 



from which one can define 



1,2,.. 



I4{'m -i— n) — U„im-iUm-lm-2 ' ' ' Un+1 n 

W =Zi(7V^ 0), 



(C7) 

(C8) 
(C9) 



which are also vxi/ matrices. Furthermore, denote {/^ } 
collectively as a z^-component column vector 



— (/fe„ I fi 



2 



(CIO) 



Periodicity in the n index is understood, viz., fg — f^r, 
UiQ — UiN- The eigenvalue problem can now be cast 
into 

Unn-ltn-l = A f „ =^ W f AT = A^Fat . (Cll) 

One then diagonalizes VV with eigenvectors {g^u}. 



W 



1,2,. 



(C12) 
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and obtain the solutions to eqn. Cll 



A„,/ = (p™)''^/'" exp + l)Sk} , (C13) 

Uin ^ 0) 



(C14) 



where the integer / = 0, 1, . . . — 1 labels the unit cell 
around which the Wannier states |$a„ /) are localized. 

In the continuum limit iV — >■ cx), the basis-overlap ma- 
trix Unn-i is related to the non-Abelian Berry connec- 
tion matrix A via 



ab 



where Aab = {ipk\ idk , hence 



U{m ^ n) ^ Pexp <( z / dk A{k) 



W Pexp { i J dkAik) 



(C15) 



(C16) 



(C17) 



where P denotes path ordering. VV is nothing but the 
Wilson loop operator, and is now unitary, i.e., pw 1- 
(X) shares the same eigenstates with GRQ, with 
its eigenvalues given by the phases in eqn. |C13| 



27r 



(C18) 



This agrees with the continuum limit resullP^. We shall 
take this as the Wannier centers even when N is finite. 

Eqn. |C8| would seem to suggest that the Wilson loop 
phases depend on the phase convention of {ji/'fe )} at all 
kn, but this is not true. Notice that 



(C19) 



where Pk^ = J2c=i IV'^^) (V'fc„ I is the internal space pro- 
jector onto occupied bands at fc„, and does not depend 
on the phase convention. Thus 

[U{m ^ n)]^^ = (^LIU(™ ^ ' (C20) 

U(m ^ n) ^ n.,„n.„_,^.„,_. • • • Pk^^.Pk^ (C21) 



U can be understood as a g x g matrix in the sublat- 
tice basis (recall that q is the total number of bands). 
Correspondingly, 



Pi Pi Pi 



p p p 



(C22) 



and VV is a X submatrix of V^WV in the occupied- 
band block at fcjv, where V diagonalizes H{ki,k2 — 2tt). 
In fact 



W 




v occupied bands 
q — v empty bands 



(C23) 



since any matrix element involving the empty bands are 
projected out by the head and tail Pfc„ . Thus W and W 
have the same non-zero eigenvalues. 

Appendix D: Invariance of eigenvalues of product of 
projectors under cyclic permutation 

Assume 



PiP2---Pn\'^)=X\'>P) , A^O 



(Dl) 



where {Pn\ is a set of arbitrary projectors (square ma- 
trices zero-padded to a common dimension if necessary) . 
Now apply the tail projector to both sides, one gets that 



PnPiP2 - ■ ■ PN-liPN^)) ^ KPNm , 



(D2) 



i.e., Pjs! is an eigenvector of the permuted product 
PnPiP2 ■ ■ ■ Pn-1 with the same eigenvalue. Carrying out 
the procedure recursively, we have that 

PnPn+l ■ ■ ■ PnPiP2 ' ' ' Pn = X^n) , (D3) 

\^P„} = P^Pr^+i ■ ■ ■ Pn W , (D4) 

i.e., Itpn) is an eigenstate with the same eigenvalue A af- 
ter cyclicly permuting the product for n times. Since all 
non-zero eigenvalues are preserved and cyclic permuta- 
tion does not change dimension, we have that all eigen- 
values are preserved. 

A corollary is that for any two projectors P and 
R, PRP and RPR have the same spectrum, because 
spec{PRP) = spec{RPP) = spec(i?i?P) = spec(i?Pi?). 
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